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The efficiency of the spiral-in of a black hole to the Galactic centre 
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ABSTRACT 

We study the efficiency at which a black hole or dense star cluster spirals in to the Galactic 
centre. This process takes place on a dynamical friction time scale, which depends on the value 
of the Coulomb logarithm (In A). We determine the accurate value of this parameter using the 
direct A^-body method, a tree algorithm and a particle-mesh technique with up to 2 million 
plus one particles. The three different techniques are in excellent agreement. Our result for the 
Coulomb logarithm appears to be independent of the number of particles. We conclude that 
In A = 6.6 ± 0.6 for a massive point particle in the inner few parsec of the Galactic bulge. 
For an extended object, like a dense star cluster, In A is smaller, with a value of the logarithm 
argument A inversely proportional to the object size. 

Key words: methods: N-body simulations - methods: numerical - celestial mechanics, stellar 
dynamics - black hole physics - Galaxy: bulge. 



1 INTRODUCTION 

The very young objects near the Galactic centre, such as th e 
Quintuplet star clu ster jNagata et alj Il990t lOkuda et all ll99(J) . 
the Arches cluster |N a£ata_etanE?95l) and the ce ntral star clus- 
ter flamblvn & Riekd 1 19931 iKrabbe eta?] 1 19951) are of con- 
siderable interest. One of the more interesti ng conundrums is 
the presence of stars as young as few Myr l Tamblvn & Riekd 
1 1993tlKrabbeetal.il 19951) within a parsec from the Galactic cen- 
tre lGerhardl200ll) . In situ formation is problematic, due to the 
strong tidal field of the Galaxy, which makes this region inhos- 
pitable f or star formation. One possible solution is provided by 
iGerharol 1200 ll) . who proposes that a 10 6 M© star cluster spirals 
in to the Galactic centre within a few million years from a dis- 
tance > 30 pc. The infall process is driven by dynamical fric- 
tion ( Chandrasekhar 1943). A quantitative analysis of this model 
by McMillan & Portegies Zwart 1 2002) confirms Gerhard's results. 
The main uncertainty in the efficiency of dynamical friction, and 
therewith the time scale for spiral-in, is hidden in a single parame- 
ter, called the Coulomb logarithm In A. Accurate determination of 
this parameter is crucial for understanding this process. Neverthe- 
less, a precise value of In A for the Galactic central region is not 
available. In this paper we determine In A for the Galactic centre. 
We focus on the efficiency of the interaction between an interme- 
diate mass black hole (BH hereafter) and the stars in the Galactic 
central region. In Section|3]we comment on how this approach can 
be applied to star clusters that sink to the Galac t ic cen tre. 

In the classical study of [chandrasekhar (1943), dynamical 
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friction is driven by the drag force experienced by a point mass 
moving through an infinite medium of homogeneous density. The 
consequences of finiteness and non -homogeneity have been anal- 
ysed i n various works (s eelMaoJ 1 19931 iMilosavlievic & Merrittl 
2001). lust & Penarrubia (2003) carried out an analytical study of 
dynamical friction in inhomogeneous systems, leading to a value 
of the Coulomb logari thm t hat de pends on the infalling object po- 
sition. IColpi & Pallavicinil ll998l) developed a general theoretical 
framework for the interaction of a satellite with a primary galaxy, 
able to describe dynamical friction in finite, inhomogeneous sys- 
tems. They applied their theory of linear respon se to orbital decay 
of satellites onto a spherical galaxy l Colni 1998) and short-lived en- 
counters with a high-speed secondary IColpi & Pallavicinilll998l) . 
They studied evolution of sate llites in isothermal spherical haloes 
with cores IColpi et alJll999h . extended in iTaffoni etal] <2003ft . 
treating satellite finite size and mass loss. Still, the original expres- 
sion of Chandrasekhar is used to model dynamical friction in many 
astronomical situations (see Binney & Tremaine 1987, Section 7.1; 
lHashimoto et aJ2002l) . The cases we study here are characterised 
by a point mass, with a very small mass compared to the primary 
system. Therefore Chandrasekhar' s formulation is appropriate in 
our cases. 

Dynamical friction is important for a large variety of astro- 
nomical phenomena, e.g. plan et migration ( Goldreich & Tremaine 
ll98CtlCionco & Brunini2 002). core collapse in dense star clusters 
( Portegies Zw^rt_etaLjl999j) or mergers in galaxy clusters (Makino 
ll997tlCora et all 1997t Ivan den Bosch et alJl999l) . The physics of 
the infall process of a satellite in the parent galaxy is basically the 
same as in the case of a BH spiralling-in to the Galactic centre. The 
relevant parameters, however, are quite different in the two cases. 
For example, an inspiraling galaxy has finite size, whereas a BH is 
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a point mass. Dynamical friction also plays an important role in the 
evolution of the black hole binary formed after the merging of two 
galax ies both hosting a BH at their centre iMilosavlievic & Merrittl 



2001). In this case, dynamical friction is important in the early 
phase of galaxy merging, when the BHs orbits converge and be- 
come bound. 

We determine the value of In A for a BH spiralling-in to the 
Galactic centre by means of self-consistent TV-body simulations. 
This is by far not an easy task. TV-body models either lack in the 
number of particles (a direct TV-body code can treat up to about 
10 5 particles, com pared to 10 8 f or the real system) or have to in- 
troduce softening (Aarseth 1963) and approximation of the force 
calculation (treecode (Barnes & Hut 1986) or particle-mesh code 
(Hock nev & Eastwoodll988l) ). The softening parameter e was in- 
troduced to limit the strength of the mutual gravitational interac- 
tion during close stellar encounters. Without softening, the very 
high accelerations experienced by the encountering bodies would 
cause very tiny integration steps, which would result in a substan- 
tial freeze of the global system evolution, with consequent dramatic 
performance degradation. The use of such approximation does not 
invalidate the numerical results, as long as the simulated system 
is studied on a time scale shorter than the relaxation time scale 
iBinnev & Tremainel 19871 ch.4, see also discussion in Section l3~3l 
below). The dynamical friction time scale of the systems we simu- 
late is in all cases shorter than the relaxation time scale, so we can 
safely use the approximate methods. 

Nevertheless, since close encounters have a relevant effect on 
dynamical friction, decreasing their strength by means of softening 
also decreases the strength of dynamical friction, i.e. lowers the 
value of In A. The same role of softening is played, in the particle- 
mesh code, by the grid cell size I. 

Our methodological approach for the present work (see Fig.Q 
consists of comparing the "exact" results obtained with the direct 
method for low particle numbers (up to 80 000) with the results of 
the treecode, which are less accurate and influenced by force soft- 
ening, to understand how the softening e influences the results and 
how they have to be scaled according to the value of e. Then the re- 
sults of the treecode are compared to the results of the particle-mesh 
code, to see how softening (tree) and grid-resolution I (particle- 
mesh) can be compared and scaled. Finally, having the right scaling 
between the different codes, we will be able to perform high parti- 
cle number simulations (up to 4 ■ 10 7 ) with the particle-mesh code 
to obtain the value of the Coulomb logarithm for the inner Galactic 
Bulge. 



2 METHODS AND MODEL 
2.1 Direct method 

For our direct TV-body calculations we used the kira 
integ rator module of the Starlab software environment 1 
IPortegies Zwart et alj 1200 ll) . Conceived and written as an 
independent al t ernativ e to Aarseth's NBODY4 and NBODY5 
jAarsethl I1985L Il999l) . the workhorses of collisional TV-body 
calculations for the past 25 years, kira is a high-order predictor- 
corrector scheme designed for simulations of collisional stellar 
systems. This integrator incorp orates a Hermite integration 
scheme ( M akino & Aarsethl 19921) and a block time step scheduler 
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Figure 1. A sketch of the strategy that we adopt in order to explore the e-N 
parameter space. 



(McMillan 1986) that allows homogeneous treatment of all objects 
in the system. 

While kira is designed to operate efficiently on general- 
purpose computers, it achieves by far its greatest speed when 
comb ined with GRAPE-6 special purpose hardware l Maki no et alJ 
2002 f. For the work presented here we performed simulations 
with the GRAPE-6 system at the University of Tokyo with up to 
80 000 particles. 



2.2 Treecode 



The hierarchical treecode ( Barnes &Hutll 19861) is widely used for 
the simulation of collisionless systems. The force on a given parti- 
cle i is computed by considering particle groups of larger and larger 
size as their distance from i increases. Force contributions from 
such groups are evaluated by using truncated multipole expansions. 
The grouping is based on a hierarchical tree data structure. Such hi- 
erarchical tree is realised by inserting the particles one by one into 
the initially empty simulation cube. Each time two particles are into 
the same cube, this is divided into eight 'child' cubes, whose linear 
size is one half of its parent's. This procedure is repeated until each 
particle finds itself into a different cube. Hierarchically connecting 
such cubic cells according to their parental relation leads to the hi- 
erarchical tree data structure. When force on particle i is computed, 
the tree is traversed looking for cells that satisfy an appropriate ac- 
ceptance crit erion; if the cell is not acc epted, then its children are 
checked. See lSalmon & Warrenl 1 19941) for a detailed overview on 
acceptance criteria. By applying this procedure recursively starting 
from the tree root, i.e. the cell containing the whole system, all the 
cells satisfying the acceptance criterion are found. 

Our treecode simulations were initially performed with both 
a cod e written by Jun M akino lMak inolll99ll) . and with GAD- 
GET ISprineel et alfcoOlh . In GADGET each particle is assigned 
an individual time-step, and at each iteration only those particles 
having an update time below a certain time are selected for force 
evaluation. This criterion was originally introduced in the direct 
TV-bo dy code (cf. | Aar seth 1999). This code is parallelized using 
MPI fMessagePa ssing Inte rface ForumllT997l) . In the parallel ver- 
sion, the geometrical domain is partitioned, and each processor 



1 See: |http : / /manybody ■ org| 



2 See: http : / /www . astrogrape ■ org| 



The efficiency of the spiral-in of a black hole to the Galactic centre 3 



simulated area (see Fig. |2j- This provides higher resolution only 
where it is necessary. Our PM-code is parallelized using MPI. 
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Figure 2. The different grids of SUPERBOX for a number of cells per di- 
mension n = 4. The finest and intermediate grids are focussed on the object 
of interest. Each grid is surrounded by a layer of two halo cells. Such haloes 
are not shown here. 



hosts the particles located in the domain partition assigned to it. 
The computation of forces on the selected i-particles is performed 
by scattering the particle data to remote processors. Then partial 
forces from the particles hosted by the remote processors are com- 
puted locally. Finally, calculated forces are received back by the 
i-particle host, and added up resulting in the total force on the i- 
particles. 



2.3 Particle-mesh code 

To perform calculations using several millions of partic les we 
use a particle-mesh (PM) code named SUPERBOX l Fellha uer et alJ 
2000). With the particle-mesh technique densities are derived on 
Cartesian grids. Using a fast Fourier transform algorithm these 
densities are converted into a grid-based potential. Forces acting 
on the particles are calculated using these grid-based potentials, 
making the code nearly collisionless. The current implementation 
completely neglects two-body relaxation causing it to retain only a 
small amount of grid-based relaxation jFellhauer et all200d) . 

The adopted implementation incorporates some differences to 
standard PM-codes. State of the art PM-codes are using a cloud- 
in-cell (CIC) scheme to assign the masses of the particles to the 
grid-cells; which means that the mass of a particle i is split up into 
neighbouring cells according to its distance to the centre of the cell. 
Forces are then calculated by adding up the same fractions of the 
forces from all cells to particle i. In contrast our code uses the "old- 
fashioned" nearest-grid-point scheme, where the total mass of the 
particle is assigned to the grid-cell the particle is located in. Forces 
acting on the particle are then calculated only from the forces acting 
on this particular cell. To achieve similar precision as CIC, we use 
space derivatives up to the second order. 

To achieve high resolution at the places of interest, we incor- 
porate for every simulated object (e.g. each galaxy and/or star clus- 
ter or disc, bulge and halo) two levels of sub-grids which stay fo- 
cused on the objects of interest while they are moving through the 



2.4 The theory of the Coulomb logarithm 

Dynamical friction affects a mass moving in a background sea 
of lower mass objects. A practical expression for the strength 
of the drag force on a point particle with mass Mbh is 
iBinnev & Tremainel 19871 p. 424): 
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Here X = vbh l(\Plo), a being the Maxwellian velocity disper- 
sion, and p is the background st ellar density. 

The classical value of A is iBinnev & Tremainel 1987L p. 423) 
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Here b max is the largest possible impact parameter for an encounter 
between the massive point particle and a member of the background 
population, v typ is the typical speed of the objects in the back- 
ground population, and m is the mass of each of the background 
stars. Equation J2j can then be generalised to 



A = 



(3) 



Here b m i n is the distance below which an encountering particle 
is captured, instead of being scattered by the massive object. It is 
somewhat smaller than the 90° turn-around distance. With the di- 
rect iV-body technique, A can be measured precisely. However, 
with approximate TV-body methods, such as the treecode or the 
PM code, we have to take care of the interference of the soften- 
ing length/cell size with b m i n , as discussed in Section l2~5l 

iMcMillan & Portegies Zwartl 120021) obtained an analytic ex- 
pression for the distance r(t) of the BH to the Galactic centre, with 
the assumptions that the BH's orbits are nearly circular, and the 
mass profile of the Galaxy is given by a power law: 



M{R) = AR a 
They obtained: 
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a bein g the velocity dispersion. In (McMillan & Portegies Zwartl 
(2002) the value of X in the Galactic centre is also computed, 
resulting in X — \J2 — a. Finally, we take Ro equal to the half- 
mass radius Rhm (see Section l2~6l . The best fit of equation J5J on 
the simulation data gives the value of In A for that simulation. Such 
values, for all simulation performed, are reported in the last column 
of Tables|3l@and|5l 



2.5 



The role of softening in the determination of the 
Coulomb logarithm 



As already mentioned in the introduction, softening was introduced 
in numerical stellar dynamics to limit the strength of mutual forces 
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during close stellar encounters, mainly for computational perfor- 
mance purposes. It consists in a modification of the Newton law 
for the gravity exerted by a particle j on a particle i, as follows: 



= G 



(r 2 - 



+ £ 2)(3/2) 



(6) 



0, 



where Vij = r_, — , and e is the softening parameter. As r; 
the presence of e causes the force to change from inverse square to 
elastic, with constant Gmiirij/e 3 . In this way the strength of the 
mutual force between encountering particles is no more unbound. 

Softening was first used by Aarseth 1 1963) in a particle- 
particle (PP) context. Accuracy requirements soon led to a more 
precise treatment o f close encounters a nd binaries by mean s of an 
analytic approac h ( Kustaanheimo & Stiefei [T965t lAarsefhlll972l 
iMikkola & Aarsethlll990l) . The softened force in equation ^6) re- 
mains used in the treecode, where high accuracy in close encoun- 
ters treatment is not mandatory. Here we will use the softening both 
in the treecode simulations, where it is necessary, and in the PP 
code simulations, where it is used to compare the results of the two 
codes, in order to study the relation between e and In A. 

For the PM code, as described in Section l2~3l force is not com- 
puted by using the Newton force, or the softened force in equa- 
tion {6). Instead, the fact that the gravitational potential on each 
grid point of the mesh is obtained from a density field defined on 
the same mesh, leads to an accuracy for the force on each particle 
limited by the cell size of the grid, /. 

Here, we are concerned with the accuracy of the computation 
of the encounters experienced by a black hole spiralling-in to the 
Galactic centre. Since the softening (PP and treecode) and the cell 
size (PM code) affect this accuracy, we will use e and I to quantify 
the accuracy decrement in our simulations. In Section l3~4l we will 
study quantitatively the dependence of In A on e and I. 

The reference value for e in the work presented here will be 
en = 0.003735 (units given belo w in Table0. This value, accord- 
ing to lAthanassoula et alj l2000l) . is of the same order of ma gnitude 
as the optimal softening for a Dehnen sphere distribution l Dehnen 
119931) . This distribution is similar to the power law distribution that 
we use, at least for what concerns the high central density peak, 
which is the key physical factor in the determination of the opti- 
mal softening. For an 80000 particle distribution, en is about 15 
times smaller than the mean inter-particle distance i at the initial 
BH position 7?o — 0.87 (see Section l2~oV This value for e is small 
enough to avoid spurious effects in the force between a star and its 
neighbours, but is sufficient to inhibit very close encounters. The 
expression for I can be obtained as: 
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where n is the star number density, and 
1 dM Aa 



-R" 



1 4nR 2 dR 4tt 
We used the expression in equation J4) for M , and the fact that the 
N stars in the system have the same mass m = 1/N. 

One of the effects of softening is a damping in the BH infall 
at very small values of the galactocentric distance, more notice- 
able as N increases. This can be explained with the fact that the 
inter-particle distance I decreases as the BH approaches the Galac- 
tic centre (see equation^. When £ becomes comparable to 2e, the 
role of softening in the force equation becomes dominant, since 
particles begin to 'overlap'. With N = 400 000, we get t — 2eo 
when R ~ 0.064, which is close to the value at which the damping 
arises, as Fie. ll2l below clearly shows. 
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Figure 3. Particle ranges for the simulations performed by each method. 
Crosses denote the particle values used. 



2.6 Initial condition 

We generate the initial mass distribution according to the power law 
given by equation J4j, with a = 1.2, which reprodu ces the mass 
distrib ution in the centre of the Galaxy, according to Mezg eret alJ 
( 1999). The scale factor is A = 4.25 • 10 6 Mq, corresponding t o 
0.44 in the A-body standard units jHeggie & Mathieul Il985h . 
which are reported in TableQ We use the standard units hereafter, 
unless other units are explicitly reported. The distributions that we 
generate are truncated at R = 1.7 = 13.6 pc, with a total mass within 
this radius Mtot = 1. The particles have equal mass m. Particles 
are assigned Maxwellian velocities, then the system is virialised to 
dynamical equilibrium. Then, before inserting the black hole (BH) 
particle, we let the system evolve for a few crossing times. The sys- 
tem reaches a stable configuration, whose mass profile is no more 
perfectly reproduced by Eq.|4| The best fit for A and a on the mass 
profile of the stable configuration gives: 



A = 0.53 
a = 0.9 . 



(8) 



In fact, the mass profile having these coefficients diverges 
from the original one as the distance R increases. On the other 
hand, in the region R < 2, where we study the BH infall, the 
discrepancy between the two mass profiles is small. The relaxed 
profile values are within 10% of the initial profile values. Nev- 
ertheless, for consistency we will use the values in equation {8j 
for A and a hereafter. This results in values of In A ~ 10% 
smaller than the ones given by a mass profile with coefficients 
a = 1.2 and A = 0.44. 

The BH particle is placed at the half-mass radius Rhm — 0.87 
with a circular orbit velocity, and its mass is Mbh = 0.000528. 
The background particles number varies from 16000 to 2 mil- 
lion. The low particle number simulations are performed with the 
PP code, the intermediate and high number simulations with the 
treecode and the PM code. Fig.|3|shows the particle range for each 
code. This allows us to span a large range in particle number, so that 
the influence of granularity in the BH motion towards the Galaxy 
centre can be studied. 

In contrast to the other models, we choose physical units for 
the PM code simulations. The conversion factors from physical 
units to A-body units are shown in Table Q where l.u. denotes 
the unit length in A-body units, m.u. the unit mass, v.u. the unit 
velocity and t.u. the unit time. 

The parameters of the PM calculations are chosen in the fol- 
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Table 1. Conversion table between the TV-body units used in our work, and 
physical units. The iV-body units are such that G = 1, Mtot = 1, and 
E tot = -0.25. 

G = 1 v.u. 2 l.u./m.u. 

= 4.3007 ■ 10~ 3 km 2 pc/s 2 M Q 
= 4.49842 ■ 10" 3 pc 3 /Myr 2 M© 



■ 1 m.u. 



Table 3. Overview of the PP runs. ./V is the number of particles, e is the 
softening parameter, eo = 0.003735, Mg^/m is the ratio between the 
BH mass and a particle mass, and e / b m i n the ratio between the softening 
parameter and the minimal impact parameter. 



N 




M BH /m 




In A 


16K 





8.5 





3.8 


16K 


1 


8.5 


2.6 


3.6 



42.3 6.6 

42.3 0.03 6.0 

42.3 0.3 5.3 

42.3 2.6 4.8 

42.3 5.3 3.5 

42.3 21.2 2.8 

42.3 42.4 1.8 



1 l.u. 
1 m.u. 

1 v.u. 
1 t.u. 



8pc 

1.18 ■ 10 s M 



G ■ 1 m.u. 



1 l.u. 



ll.u.' 



251.86 km/s 



■ 0.031 Myr 



80K 
80K 
80K 
80K 
80K 
80K 
80K 





0.01 
0.1 

1 

2 
8 

16 



Table 2. Resolutions (i.e. cell sizes) of the different grid levels for the dif- 
ferent choices of n in the PM code, n denotes the number of cells per di- 
mension. The cell sizes of the different grid-levels (outer, middle and inner) 
are given in pc. 



n 


outer 


middle 


inner 


32 


10.00 


0.69 


0.17 


61 


4.67 


0.32 


0.08 


128 


2.26 


0.15 


0.04 



lowing way: the grid sizes are kept constant at 

-Rsystcm = 140.0 PC 

Rout = 9.6 pc (9) 

-Rcoro = 2.4 pC 

and are focussed on the center of mass of the "bulge" model, as 
sketched in Fig.|2| To change the resolution we alter the number of 
grid cells per dimension from 32 up to 128. With this choice the 
cell sizes listed in Tableware achieved. 

To avoid a self-acceleration of the black hole we choose 
Ho,o,o = 1-0 of the Greens-function, which worse ns the energy 
conservation of the code (see iFellhauer et al]|20 00) especially in 
simulations with high particle numbers per grid cell, but suppresses 
the forces of particles in the same cell on each other and on them- 
selves. 

To speed up the simulations, the time step in the PM code sim- 
ulations should be as large as possible, but small enough to prevent 
spurious results. Therefore we started with a time step of 1,000 yr 
and reduced it to 200 and 50 yr. The results of the 200 yr and 50 yr 
time step do not differ from each other, therefore the global time 
step is chosen to be 200 yr. Conversely, the time step in the treecode 
and direct code simulations is variable and different for each parti- 
cle. Time step values are in this case in the range 2-30 000 yr, with 
about 90% of them in the range 100-300 yr. 



3 RESULTS 

We will now study the dependence of our results on the number of 
particles iV in Section l3~T1 and compare the various iV-body meth- 
ods with identical initial realisations in Section l3~2l After having 
convinced ourselves that the various techniques produce consistent 
results, we continue by studying the effect of softening/cell size 



Table 4. Overview of the treecode runs. Meaning of symbols is the same as 
in Tablel3labove. 



N 




M BH /m 




In A 


80K 


1 


42.3 


2.6 


4.7 


400K 


1 


211.3 


2.6 


5.0 


2M 


1 


1056.5 


2.6 


4.9 


80K 


0.1 


42.3 


0.3 


5.7 


80K 


2 


42.3 


5.3 


4.1 


80K 


8 


42.3 


21.2 


3.0 


80K 


16 


42.3 


42.4 


2.0 


80K 


32 


42.3 


84.7 


1.6 


80K 


1 


84.5 


1.3 


5.4 


80K 


1 


169.0 


0.7 


4.6 


400K 


1 


422.6 


1.3 


4.6 


400K 


1 


845.2 


0.7 


4.2 



Table 5. Overview of the PM runs. N is the number of particles, n the 
number of grid cells per dimension, m the particle mass, / the intermediate 
grid cell size, N c the average number of particles per cell, m c the average 
mass of a cell, Mgu / m c the ratio between the BH mass and the cell mass, 
and finally I / b m i n the ratio between the cell size and the minimal impact 
parameter. 



N 


n 


rn 
[M ] 


I 

[pc] 


l cell J 


m c 

L cell 1 


Mbh 


; 


In A 


80K 


16 


1475 


1.60 


46.3 


68287.0 


0.9 


114.3 


n/a 


80K 
400K 
2M 


32 


1475 
295 
59 


0.69 
0.69 
0.69 


3.6 
18.2 
91.1 


5375.4 
5375.4 
5375.4 
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Figure 4. Time evolution of the radial distance of the black hole to the 
Galactic centre. The various curves (identified in the top right corner) 
present results obtained with the treecode. the X-axis is presented in TV- 
body time units: one TV-body time unit corresponds to about 0.031 Myr. The 
distance of the black hole to the Galactic centre (Y-axis) is given in terms 
of its initial distance. In these simulations is e = 0.003735 ~ 0.03 pc and 
M B h = 0.000528. 

(Section f3.3t and black hole mass (Section [3.5t on the value of the 
Coulomb logarithm in the inner part of the Galaxy. 

Our simulations aimed at several goals. 1) understanding the 
scaling of the system dynamics with respect to the number of parti- 
cles N, and within this scaling, how results from different methods 
compare with each other. 2) How, at a fixed value of N, the soften- 
ing parameter influences the dynamics, changing the value of In A. 
The particle-mesh method does not make use of softening. The cell 
size in the PM code can be seen in this context as a softening length. 
In our framework, it is crucial to understand the relation between 
the PP code and treecode softening parameter and the PM code cell 
size. 3) We also study how the BH mass influences the infall time. 
We doubled and quadrupled the BH mass, and observed how this 
affects the value of In A. 

A resume of all the runs that we performed is reported in Ta- 
ble[3]for the PP code runs, Table0for the treecode runs, and finally 
Table[5]for the PM code runs. In all of our runs, the system remains 
in equilibrium during the whole BH infall, with no significant mass 
loss from stellar escapes, and a mass profile independent of time. 



3.1 Dependence of In A on TV 

In order to obtain a precise measure of In A, ideally one would run 
a direct TV-body simulation with TV of the order of the number of 
stars in the Galactic bulge, which amounts to ~ 10 8 . Such high 
number makes a direct simulation unfeasible, and imposes the use 
of approximate methods instead. In order to evaluate the reliability 
of the approximate methods, we compared the PP code runs with 
the treecode runs. The PP code runs give a reliable picture of the 
system dynamics at low particle numbers, i.e. at high granularity. 
Using the treecode we can reach a much highernumber of particles, 
up to two million, which still is two orders of magnitude lower than 
the real system. A comparison of the results from the two methods 
allows us to estimate the validity of the treecode runs, up to 2 mil- 
lion particles. Then we can compare the treecode runs and the PM 
runs, in order to validate the results from the latter, which has the 
capability to simulate systems of about 100 million stars. In such a 
way we will eventually be able to study the infall of a BH into the 
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Figure 5. Same as Fig.|4]above, but for PM code simulations. The inter- 
mediate grid cell size is here I = 0.69 pc, and Mgg = 0.000528. 
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Figure 6. Black hole infall at various cell sizes, and large difference in 
TV. Results here are from PM code simulations. The case TV = 80 000, 
I ~ 5en is not shown for readability reasons, since it would overlap with 
the I ~ 10en results. 

Galactic centre in a simulation environment with a realistic value 
of TV. 

In Fig. [4] we show the evolution of the BH distance from 
the centre of mass of the system for three treecode simula- 
tions. TV varies from 80000 to 400000 and 2 million, with 
e = eo = 0.003735, corresponding to about 0.03 pc. In Fig. [5] 
we present a similar figure from PM code simulations. Here is 
TV G {80 000,400 000,2 000 000}, with 32 cells per dimension, 
resulting in a cell size of about 0.69 pc. 

Figs |4] and [5] show that increasing TV results in a much 
smoother motion of the BH in its infall towards the centre of the 
Galaxy. The BH infall rate is not much affected by a change in TV. 
Accordingly, the value of In A for each of the two sets above is 
consistent, as values in Table[4](first three rows) and Table[5](rows 
with I = 0.69) show. 

In order to study further the extent of the influence of TV on the 
infall rate of the BH, and hence in In A, we compare in Fig.|6|results 
from PM code simulations with increasing grid refinement, and 
extreme difference in TV. To quantify the grid resolution, we use 
the cell length at intermediate refinement, which is the cell length 
pertaining to the physical region where the BH evolves for most 
of its infall. We measure this length in units of eo = 0.003735, 
which makes the comparison with the softening parameter of the 
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Figure 7. Top panel shows a black hole infall simulated by the PP code, 
with N = 80 000, M B h = 0.000528 and e = 8e . Bottom panel shows 
a comparison of the PP results with treecode and PM results. Parameter val- 
ues are in all cases the same, except for the PM cell size, which is I = 10eo. 
Plotted values are averages over 10 time units. 

treecode easier. N has no remarkable influence on the infall rate, 
except for the case where the cell size is I = 0.15 pc ~ 5eo- In 
this case the simulation with N — 80 000 ((data not reported in 
the figure), shows an incorrect BH infall, comparable to the case 
I — 0.32 ~ 10eo • This can be explained by the fact that in the low I, 
low N case, the cells are so small, and the particles so few, that 
many cells in the PM grid are empty (see also the N c column in ta- 
ble|5| which gives the average number of particles per cell). When 
N c <C 1, the density field is incorrect, with many grid points hav- 
ing a null value, because the corresponding cell is empty. In this 
condition, the gravity field computed by the PM code becomes un- 
reliable, affecting the numerical results, as in the simulation with 
N = 80 000 and I ~ 5e . 

3.2 Comparison of the codes 

In this section we compare the results obtained from the various 
codes, to check their consistency. The comparison of the PM re- 
sults with the two other codes results is particularly critical, since 
the PM code computes forces using a different mathematical ap- 
proach, i.e. a grid based force derivation vs a direct particle-particle 
computation for the PP code, or particle-multipole computation for 
the treecode. A consequence of this is a different parameter to tune 
the accuracy of the simulation, namely the cell size I for the PM 
code, and the softening length e for the other two codes. We will 
study here how these two parameters influence the black hole infall. 

In Fig.Qwe show the time evolution of the galactocentric BH 
distance R simulated by the PP code, accompanied by a plot of 
the time evolution of AR/Rpp for treecode and PM simulations, 
where AR = (R — Rpp). The relative difference AR/Rpp re- 
mains small for a large fraction of the infall, and the final discrep- 
ancy is mostly due to the small values of the quantities at that point, 
which are likely to amplify relative differences. As the following 
figures also show, the BH infall is predicted with very good consis- 
tency among the codes. 

In Fig.|8|selected treecode runs with N = 80 000 and increas- 
ing e are compared with the direct code runs having the same values 
of N and e. At the same time, the figure shows how the infall time 
increases (and implicitly how In A decreases), as e increases. Fig. [8] 
and Table |6| show that the results from the treecode, the PP code, 




Figure 8. Comparison of results from the PP code with results from the 
treecode, at different values of e. For all cases shown here is N = 80 000 
and Mbh = 0.000528. The PP simulation with e = 8eo has been already 
shown in Fig. 171 
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Figure 9. Comparison of PM results with treecode results. PM simulations 
have cell size / equal to resp. 10eo and 23eo; softening parameters in the 
treecode runs are resp. 8eo an d 16erj. In all the above cases, is N = 80 000 
and M B h = 0.000528. 



and the PM code are in good agreement. The agreement of the re- 
sults from the three methods, and the scaling of In A with e, will be 
further studied quantitatively in Section l3l4l 

In order to understand how the cell length I of the PM code and 
the softening parameter e of the PP code and treecode relate with 
each other, we compare in Fig.|5|the results from the PM code and 
treecode simulations with 80 000 particles. The BH infall as shown 
in Fig.|5|depends on the value of I or e. Remarkably, I and e seem to 
play the same role not only qualitatively, but also quantitatively: in 
a PM run, a given value of I induces an infall which is very similar 
to the infall, in a treecode run, with e assuming that same value. In 
Section l3l4l this relation will be studied further. 



3.3 The effect of softening/grid 

The influence of the softening parameter on the BH dynam- 
ics has been studied by running a number of simulations with 
the three codes. In Table |6| we report the value of In A ob- 
tained from our simulations. For the PP code and treecode sim- 
ulations, we increase e from to 32eo = 0.1195 ~ 0.96 pc. For 
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Table 6. In A versus e from PP code, treecode, and PM code runs. 
For the PP code and treecode runs is N = 80 000. For the PM code 
runs is N = 2 million. The reference value for the accuracy parameter is 
e Q = 0.003735. 



i/erj PM code 



e/eo 


PP code 


treecode 





6.6 




0.01 
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the PM code, we increase / from 2.5 eo to 23 eo- In all cases is 
Mbh = 0.528 • 10~ 3 ~ 62 3OOM . 

For the PP code and the treecode, we selected TV — 80 000 
as a suitable value. The relaxation time is for this value of TV 
t x ~ 0.1N/ In TV ■ R/vtyp — 2000, about one order of magnitude 
larger than the typical BH infall time, so that the system is collision- 
less, and we can confidently use the treecode to simulate it. With 
this choice for TV, the BH mass is Mbh / in — 42.3, (see Table[4}- 
As a cross-check, we ran two PP runs with TV = 16 000 which, as 
expected, gave incorrect results (see Table[3j. This is due to both a 
too small Mbh /m ratio, and a too short relaxation time (t x ~ 400 
in this case). We did not increase e above ~ 0.12, since at this point 
e is already much bigger than b m in (see Table|4j, and the infall time 
is now close to t x . 

For the PM code simulations, we used TV = 2 million in order 
to have enough particles to fill all the cells, even for the simula- 
tions with a small I. As Table [5] shows, for I — 0.076pc ~ 2.5eo 
the average number of particles per cell is already TV C = 0.1. Since 
a PM simulation gives incorrect results for TV C <C 1 (see also the 
discussion at the end of Section l3"TV we did not decrease / below 
2.5 e . 

The decrease of the value of In A as e or / increases is clear 
from Table |S| In the next sub-section we focus on the relation be- 
tween A and e, and provide a fitting formula for lnA(e). We use 
hereafter e to refer either to the softening length of the PP and 
treecode, or the cell size of the PM code. As shown on Fig. [5] 
and discussed above, these two parameters play the same role even 
quantitatively in affecting In A. In this respect, we refer to e as a 
generic accuracy parameter. 



3.4 Determination of In A 

We will study here the relation between e and In A. As just said 
above, in this context e will be used as the accuracy parameter, 
and it will refer to either the softening length used in the PP and 
treecode, or to the cell size in the PM code. 

A mathematical expression for the relation between e and In A 
can be found by considering how softening affects two body scat- 
tering. The role of e is to prevent too close stellar encounters. In this 
respect, the effect of introducing a softening length is to increase 
the minimal impact parameter. Hence, we can define an effective 
impact parameter b e // = b m in + e, and we modify equation to 
become: 

bmax , bmax 



In A = In- 



°eff 



In- 



(10) 
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Figure 10. In A vs e, and best fit for In A = K — ln(a + e). Values for K 
and a are given in TablelTl The inset in the figure is a magnification of the 
low e region. In all cases is Mbh = 0.000528. For the PP and treecode 
runs is N = 80 000, for the PM code runs is N = 2 000 000. Error bars 
are omitted from the PP values to improve readability. For the same reason, 
In A values for e/eo < 1 are shown only in the inset. 



Table 7. Best values for the parameters K and a, and error on In A for the 
fit of In A = K - ln(a + e). 



PP code 



treecode 



PM code 



K -0.94 ± 0.21 
a -10~ 3 0.80 ±0.28 
A (In A) 0.6 



-0.64 ±0.10 -0.59 ±0.05 
0.88 ± 0.20 0.74 ± 0.08 
0.3 0.2 



We will now fit this equation with the values reported in Table|6| In 
order to perform the fit, we change equation i 1 01 in a more suitable 
form, as follows: 



In A = In b„ 



ln(b m j n + e) = K — ln(a + e) . 



(11) 



We will refer to bmax and b m in as the theoretical values of the 
maximal and minimal impact parameters, as they can be obtained 
from equations J2j and J3), and K and a as the corresponding ex- 
perimental values obtained with the fit. 

The best fits for K and a with respect to simulation values are 
reported in Table for all codes. Such fits have been performed 
with a fixed value for Ro, i.e. the Ro = Rhm- In fact, the not per- 
fectly circular orbit of the BH results in an oscillatory behaviour 
for the BH galactocentric radius. In this case, having Ro fixed could 
not be an appropriate choice for the fit. We checked whether having 
Ro as a free parameter in the fit leads to different results in In A. 
We obtained values for In A within the error bars in Fig. 1101 and 
values for Ro within Rhm ± 0.05. We can conclude that, although 
the galactocentric BH radius does not decreases smoothly, but in 
an oscillatory fashion, having Ro fixed to the actual initial BH ra- 
dius in the simulations leads to correct fits for the value of In A. 
With respect to the PM values, a further peculiarity is that when 
the BH enters the finest grid area, i.e. approximately at R = 0.3, 
the value of / decreases (see Section l2~3l and Fig.|2j. This causes 
b e ff to become smaller, increasing the value of In A. In fact, a fit 
of the PM data limited to values of R > 0.3 gives values of In A 
systematically higher by ~ 0.3 ~ 2A(ln A). 

From the PP code value of K in Table Q we obtain for bmax 
the experimental value b max — e K ~ 0.39. This value is quite 
smaller than what one would expect. Since b m ax has the mean- 
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Figure 11. Compai'ison of In A vs e at different initial galactocentric BH 
radii. The smaller values of In A for the cases with Rq = R qm indicate 
that b m ax is influenced by the galactocentric BH radius. 

ing of the maximal impact parameter, a natural choice is to assign 
it a value of the order of the system size, which in our case would 
result in b max = 2. The maximal radius for dynamical friction in 
our system is then about one quarter of what it is customarily as- 
sumed. PM and treecode values are slightly higher, but still quite 
smaller than b max = 2. Also a is smaller than the theoretical value 
bmin = G ■ {Mbh + m)/v% p = 1.41 ■ 1CT 3 , by a factor 3. The 
a value for all codes is perfectly consistent. 

An explanation for the discrepancy between the values of 
bmax and b max is that the BH, while moving to the Galactic cen- 
tre, is off-centre with respect to the density peak (in fact the BH 
is spiralising towards it). With respect to the BH position, the den- 
sity distribution is then asymmetric. This density peak has clearly a 
greater influence on the BH dynamics, contributing more than the 
other regions of the system to the dynamical friction on the BH. 
This leads to a value of b max affected by the galactocentric BH ra- 
dius. This approach is studied in detail by Hashi moto et alJ l20Q2). 
who propose the galactocentric radius as a value for bmax in the 
context of the spiralisation of satellite galaxies. 

In our simulations, the galactocentric radius varies from 
R ~ 0.9 at the beginning of a simulation, to R ~ at the end of it. 
The value of b max that we find is within this range, and it can be 
interpreted as an order estimate of a maximal impact parameter 
that depends on the galactocentric BH radius. 

In order to explore this aspect further, we simulated the infall 
of the same BH, starting at the quarter mass radius R qm ~ 0.43, 
for e ranging from to 16erj. What we expect is a smaller value of 
t K hence smaller values of In A. All simulations are performed 
with the treecode, except for the e = case, which is simulated 
with the PP code. Our results are in Fig. ll II We can see there how 
smaller are the values of In A for the cases when the BH starts at 
the quarter mass radius. A fit on these data gives A' ~ —1.1, which 
implies b max ~ 0.33, which is smaller than the value of b max ob- 
tained for the BH starti ng from the half mass ra dius. Our findings 
support the argument of Hashimoto et al. 1 2002). 

3.5 Varying black hole mass 

We also studied the effect of a variable BH mass on the value 
of In A. We simulated, using the treecode, the infall of a BH 
of mass two times and four times larger than the default mass 
M = 0.528 • 10" 3 ~ 0.62 • 1O 5 M . We studied this infall in 
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Figure 12. Black hole infall for different values of the BH mass, and differ- 
ent values of N. Simulations are performed with the treecode. Simulation 
results are compared with the analytic solution, equation |5J, with In A ob- 
tained from equation 1111 and Table 171 The error bars at the bottom of the 
analytic curves correspond to a variance A(ln A) = ±0.3. 



both the 80 000 particles configuration, and the 400 000 particles 
configuration. In all cases, we used our standard value for e, i.e. 
en = 0.003735. In Fie. ll2l the distance r of the BH from the centre 
of mass of the system is shown for all the cases mentioned above, 
together with the Mbh = Mo cases. From equation i 1 1 i and Ta- 
bleQ the appropriate value for In A in the above cases is: 

In A = K - ln(a + e ) ± A = 

-0.64 - ln(0.00088 + 0.003735) ± 0.3 ~ 4.7 ± 0.3 . 

We also show in Fig. ll2l the analytic curve r(t), as given by equa- 
tion with In A = 4.7. An error bar gives, for each analytical 
curve, the spread corresponding to a variance A(ln A) = 0.3. 

The results shown in Fig. ^| are consistent with the hypoth- 
esis that a variation in the BH mass has a little effect in the value 
of In A. In fact, In A shows a logarithmic dependence on Mbh 
through the parameter b m in, which depends linearly on Mbh (see 
equations J3J and J3J). Assuming that also the experimental value 
a depends linearly on Mbh, we obtain In A ~ 4.6 ±0.3, and 
In A ~ 4.4 ± 0.3 respectively for the 2 Mbh case and the 4 Mbh 
case. This results in a small displacement towards the right of the 
corresponding analytic curves in Fig. 1121 which does not affect 
the conclusions that can be drawn from the figure. The theoreti- 
cal curve fits very well with the M = 2 Mo, N = 400 000 case. 
The other simulation curves are within, or very close to, the error 
in r(t) associated to the error in In A. We can conclude that a vari- 
ation in the mass of the infalling object has little influence in the 
value of In A, which is important in view of extending this work to 
the case of the infall of a star cluster. 

The fitting formula for In A vs e was obtained from simula- 
tions with Mbh = Mo- This formula predicts In A for the cases 
with Mbh > Mo with a very good accuracy, showing that it can 
be applied in a more general context, in order to forecast the value 
of the Coulomb logarithm. 

Fig. 1121 also shows a damping in the BH infall at very small 
values of R, especially for the N = 400 000 case. This effect, de- 
scribed in Section l2"31 is clearer in the iV = 400 000 case, since the 
particle density is higher in this case, compared to the A^ = 80 000 
case. 
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3.6 Comparison with related work 

iMilosavlievic & Merrittl 1200 ll) study the dynamical evolution of 
two black holes, each one at the centre of a power law cusped 
galaxy core. They simulate the merging of the two galaxies, which 
results in the two black holes forming a hard binary at the centre 
of the merged galaxy. In Section 3 of their paper they discuss the 
value of In A in their simulations. They measure the decay rate of 
the two black holes, and compare this value with theoretical es- 
timates. When they compare their experimental decay rate with an 
estimate for the case of the infall of an isolated black hole, they find 
a theoretical estimate about 6 times lower than the measured value, 
under the assumption that In A ~ 1.6. If the value of In A is not 
theoretically pre-determined, and is instead obtained from the de- 
cay rate equation, the result is In A ~ 10. Similarly, they compare 
the experimental value with an estimate for the case of two mu- 
tually spiralling spherical distributions of matter. In this case they 
assume In A ~ 1.0, and obtain an estimate for the decay rate about 
a factor of 2 lower than the observed value. Determining In A from 
the measurement would give in this case In A ~ 1.87. The values 
of In A that we find are between the two values above. 

The value for In A ~ 1 that they assume in their theoretical es- 
timates, comes from a derivation that they present in appendix A of 
the same work. This derivation is based on results of Maoz 1 1993). 
Under the assumption that the stellar density obeys a power law 
centered on the BH position: 

p{r) = Po (j^—) , (12) 

they obtain A ~ 1/a ~ 1, which actually implies b max — b m in, 
whereas it is customary to consider b max 2> & m j ra . 

Their assumption in fact is valid only when the BH is close 
to the centre of the power law distribution. In their context this is 
true when: 1) the separation between the two BHs is much larger 
than the half mass radius of the two galaxies. In this case each BH 
is at the centre of its own galaxy, and at the same time its motion is 
not yet heavily perturbed by the other galaxy. 2) the BH binary has 
hardened, and occupies the centre of the merged galaxy. 

During the transient phase, when the two BHs have not yet 
formed a binary, the density distribution that affects the motion of 
the BHs is double-cusped, with a BH in each of the two cusps. This 
is substantially different from the density distribution modelled by 
equation 1 1 21 . 

This qualitative argument would make the density distribution 
in equation j 1 21 inapplicable during the transient phase, and could 
explain why Milosavlievic & Merritt 1 2001) find a higher than ex- 
pected value of In A in the transient. The analytical evaluation of 
In A according to the technique used by them is by no means triv- 
ial, when symmetry arguments cannot be straightforwardly applied. 
We will address this issue in future developme nts of the present 
work; the theory of linear response of IColp i & Pallavicini 1 1998) 
could be very useful in this context. 



l2002t) . lMorrisl 1 19931) proposed that a star cluster at some distance 
from the Galactic centre could spiral-in due to dynamical friction 
(see also lGerhardl <200ll) ). The efficiency of dynamical friction de- 
pends sensitively on the actual value of the Coulomb logarithm 
In A. 

4.1 Sinking of massive black holes in the Galactic centre 

We performed A^-body simulations for a large range of conditions. 
In Section l3~T1 we varied the number of particles, in Section l3~3l we 
varied the size of the object, and in Section l3~5l we varied its mass. 
With direct A^-body simulations we measured the actual value of 
the Coulomb logarithm In A. We study the behaviour of In A for 
various types of A^-body solvers and particle numbers. We also 
study the behaviour of In A as a function of the softening length 
e. Only the direct A^-body code can perform a true measurement 
of the Coulomb logarithm, because it is able to resolve even the 
smallest length scales and time scales. This, however, makes the di- 
rect code very slow and, even using the very fast GRAPE-6 special 
purpose device, we are able to perform simulations with only 10 5 
particles. This is a small number compared to the actual number of 
stars in the Galactic centre. With approximate methods (treecode 
and particle-mesh) we are able to increase the number of particles 
up to 2 million. The cost of this is a lower accuracy in calculating 
stellar motion below a typical length scale e. We studied how this 
length scale influences In A, by affecting the value of the minimal 
impact parameter. 

4.2 Young dense clusters in the Galactic centre 

The study of the dependence of In A on e described above is also 
of astronomical interest, because e can be interpreted as the typical 
length of a finite size infalling object. Based on this, our analysis of 
the dependence of In A on e can be seen as a first approach to the 
study of the infall of a star cluster of typical size e toward the Galac- 
tic centre. We found (see Fig. 1 101 that the value of In A decreases 
quite rapidly as e increases, with the logarithm argument A oc 1/e. 
The typical size of the compact young cl usters observed in the 
Galactic bulge is ~ 0.3 pc iFiger et alll 9991. which corresponds to 
e ~ 10eo- With this value of e, from equation 1 1 1 i and TableQ we 
obtain In A ~ 2.9, about 60% less than the value for a point mass. 
The infall time is roughly doubled. For our choice of object mass, 
M ~ 62 300 Mq, and initial galactocentric radius, Rq ~ 7 pc, we 
have an infall time that increases from ~ 6 Myr for the point mass, 
to ~ 12.5 Myr for an object of typical size ~ 10eo — 0.3 pc. 

We also studied the uncertainty associated with the maximal 
impact parameter b max . We found that for an infall to the Galac- 
tic centre, the infalling object is mostly influenced by the density 
peak at the Galactic centre itself. A good choice for b max is then 
bmax — (3Ro, where Ro is the initial galactocentric radius, and 
/3~0.5 



4 APPLICATIONS TO STAR CLUSTERS 

Recent observations of the Galactic Centre have revealed a pop- 
ulation of very young clusters with ages less than 10 Myr. The 
presence of such stars inside the inner pc of the Galaxy is puz- 
zling, as the strong tidal field in the Galactic centre easily prevents 
star formation. The origin of these stars is therefore vividly debated 
jGerharj200ll:lMcMillan & Portegies Zwarl2002tlKim & Morris! 



5 CONCLUSIONS 

We simulated the evolution of a massive particle in a sea of lighter 
particles in a self gravitating system. The main goal of this work is 
to obtain an accurate value of the Coulomb logarithm (In A). This 
helps us to understand the dynamics of the Galactic bulge and the 
rate at which intermediate mass black holes sink to the Galactic 
centre. We also study the effect of the finite size of the inspiraling 
object. 
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We ran both A-body particle-particle (PP) simulations, soft- 
ened treecode simulations, and particle-mesh (PM) simulations. 
The comparative simulations are performed for 80000 particles, 
and all result in the same value of In A. For a point particle near the 
Galactic centre we find In A = 6.6 ± 0.6. In addition we measure 
the change in the Coulomb logarithm with respect to the softening 
parameter e, which reveals A oc 1/e. We argue that e can be inter- 
preted as the typical length of a finite size object, such as a star 
cluster, so that In A as a function of e can be seen as a first approx- 
imation of the dependence of the Coulomb logarithm on the size of 
an infalling star cluster. 

We also observed a value of the maximal impact parameter 
bmax different from the customarily assumed value, which is pro- 
portional to the system size. We found that our results are more 
consistent to a value of b max linearly dependent on the BH galacto- 
centric radius, which is in agreement with Hashimoto et alJf2002l) . 

We performed simulations with up to 2 million particles us- 
ing a treecode. The obtained value of In A does not depend on 
the number of particles. Apparently, 80000 particles is already 
enough to eliminate any granularity for our choice of initial con- 
ditions. The results of the treecode, at the low A-limit, are in 
excellent agreement with the PP simulations, and we find the 
same scaling with respect to e. Increasing the black hole mass re- 
duces the time scale for spi ral-in as was expected from theory (see 
McMillan & Po rtegies Zwartl2002h . 

Finally we compared the results of our PP and treecode simu- 
lations with a particle mesh (PM) method. We compared the meth- 
ods for N up to 2 million. The results of our PP, treecode, and PM 
calculations are in good agreement. The cell size in the PM model 
is directly comparable to the softening length e in the PP and tree 
methods. 
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